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This  work  was  prepared  by  the  Lockheed  Palo  Alto  Research 
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F33615-69-C-1523.  It  was  administered  under  the  Structures  Division, 
Air  Force  Flight  Dynamics  Laboratory,  with  Mr.  T.  N.  3erne>in  (FBR) 
acting  as  Project  Engineer. 

This  report  was  completed  in  January  1972  and  covers  work  per¬ 
formed  between  April  1969  and  January  1972.  The  supervision  of  this 
project  was  carried  out  by  Mr.  B.  O.  Almroth  of  the  Structural 
Mechanics  Laboratory,  LMSC. 
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•^This  report  presents  a  theory  for  nonlinear  collapse  analysis  of  shells  with  gene-al 
shape.  The  theory  combines  energy  principals  and  finite  difference  methods  to 
obtain  a  system  of  nonlinear  equations;  these  are  solved  by  a  modified  Newton- 
Raphson  techrique.  For  greater  economy  and  flexibility  in  the  analysis  a 
capability  is  provided  for  use  of  variable  spacing  finite  difference  grids.  Inelastic 
material  behavior,  as  predicted  by  *he  White-Besseling  Theory,  is  incorporated 
into  the  analysis.  A  computer  code,  STAGS,  based  on  the  theory  has  been  written 
and  used  to  solve  a  number  of  sample  problems.  Results  for  these  problems 
are  presented. 
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Section  1.  0 
INTRODUCTION 


Recent  improvements  in  computer  technology  and  numerical 
analysis  methods  have  led  to  significant  advances  in  structural  analysis 
capability.  Computer  programs  are  now  available  for  analysis  of  the 
static  behavior  (linear  or  nonlinear)  of  almost  any  shell  of  revolution  pub-* 
jected  to  axisymmetric  loading.  For  nonsymmetrical  loading  or  for  shells 
of  general  shape,  a  static  analysis  is  readily  performed  provided  that  the 
response  is  linear.  This  capability  is  substantiated  in  Reference  1.  How¬ 
ever,  nonlinear  effects  are  frequently  important  in  shells.  Because  t.  ese 
structures  are  thin,  collapse  or  loss  of  stability  is  generally  the  critical 
mode  of  failure.  Thicker  shells  are  often  subjected  to  loads  of  such  mag¬ 
nitude  that  material  nonlinearities  become  important.  Reliable  and  depend¬ 
able  computational  systems  for  this  important  class  of  problems  have  not 
been  developed,  although  there  are  computer  codes  available  for  some 
special,  cases. 

Several  years  ago,  a  research  program  was  initiated  at  Lockheed 
with  the  goal  of  developing  a  computational  system  for  such  nonlinear 
problems.  This  program  has  resulted  in  the  STAGS  (STructural  Analysis 
of  General  Shells)  computer  code  for  analysis  of  the  static  nonlinear  re¬ 
sponse  of  general  shells.  STAGS  is  based  on  a  theory  in  which  the  shell 
surface  is  subdivided,  by  means  of  a  finite  difference  grid  work,  into  a 
set  of  subareas.  The  strain  energy  density  for  each  subarea  is  then  ex¬ 
pressed  in  terms  of  displacement  components  and  their  derivatives.  After 
the  derivatives  have  been  replaced  by  their  finite  difference  equivalents, 
the  energy  can  be  calculated  and,  together  with  the  potential  energy  due 
to  applied  loads,  summed  over  the  shell  surface.  The  total  potential 
energy  expression  of  the  shell  so  obtained  is  then  minimized  according  to 
familiar  energy  principles  and  a  system  of  nonlinear  algebraic  equations 
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in  the  unknown  displacements  results.  These  equations  are  solved  by  a 
Newton-Raphson  technique. 

STAGS  is  an  outgrowth  of  work  on  the  buckling  of  cylindrical 
panels  with  nonuniform  membrane  stresses  that  was  initiated  at  LMSC 
in  1963  under  the  sponsorship  of  NASA  Marshall  Space  Flight  Center 
(Ref.  2).  The  basic  nonlinear  computer  program  for  cylindrical  shells 
with  cutouts  (Ref.  3)  and  a  linear  version  including  analysis  of  free  vibra¬ 
tions  (Ref.  4)  were  developed  under  the  LMSC  Independent  Research  Pro¬ 
gram.  Under  contract  with  the  Naval  Ship  Research  and  Development 
Center  (NSRDC),  the  linear  version  of  the  code  was  developed  to  include 
shells  of  revolution  with  smooth  but  otherwise  arbitrary  cutouts  (Ref.  5). 

The  work  reported  here  extends  the  nonlinear  version  to  shells  of 
more  general  shapes  with  cutouts  of  arbitrary  contour.  In  addition,  inelas¬ 
tic  deformations  and  a  capability  to  handle  a  finite  difference  grid  with 
variable  nodal  point  spacing  have  been  added.  In  a  parallel  effort  funded 
by  Lockheed's  Independent  Research  Program,  the  equations  were  further 
generalized  to  include  nonorthogonal  coordinates  (Ref.  6).  As  this  work 
was  completed  before  the  end  of  the  contract  period,  it  was  possible  to 
include  the  more  general  equations  in  this  report. 

Further  expansion  of  the  STAGS  program  has  been  accomplished 
under  parallel  research  studies  funded  by  the  Air  Force  Space  and  Missile 
Systems  Organization  (SAMSO)  and  by  the  NASA  Langley  Research  Center. 
During  the  now  completed  SAMSO  study,  provisions  were  made  in  the 
STAGS  code  to  allow  both  the  temperature  and  material  properties  to  vary 
over  the  surface  and  through  the  thickness  of  the  shell.  In  addition,  a 
bifurcation  buckling  branch  was  added.;  Parameter  studies  were  made  to 
evaluate  the  applicability  of  the  bifurcation  buckling  approach  to  reentry 
vehicle  analysis  (Ref.  7).  Although  most  of  u.cse  extensions  were  made 
primarily  to  render  STAGS  suitable  for  reentry  vehicle  analysis,  they  have 
considerably  enhanced  the  overall  capability  of  the  code.. 

The  NASA  study  is  currently  in  irogre^s.  Under  this  study,  STAGS 
is  being  developed  to  handle  Begov  atei  -mi  branched  shells,  and  to  treat 


realistic  types  of  shell  wall  constructions  including  those  which  involve 
anisotropic  materials.  Finite  difference  expressions  based  on  non- 
rectangular  grids  and  an  automatic  grid  generator  are  also  being  added. 

A  time  integration  scheme  will  be  developed  and  included  in  STAGS. 

This  will  permit  the  solution  of  dynamic  response  and  dynamic  buckling 
problems.  The  NASA  work  is  scheduled  for  completion  by  the  summer  of 
1972. 

A  STAGS  user's  manual  that  documents  all  of  the  modifications 
completed  to  date  has  been  prepared  (Ref.  8). 


Section  2.  0 
STAGS  THEORY 


In  the  application  of  finite  difference  techniques  to  shell  analysis, 
it  has  been  customary  to  assume  that  lines  of  curvature  constitute  the 
surface  coordinate  lines  which  form  the  finite  difference  mesh.  This 
assumption  results  in  orthogonal  coordinate  iines  and  leads  to  simple 
shell  equations;  however,  there  is  a  serious  disadvantage  to  this  approach 
in  that  in  many  instances  shell  boundaries  do  not  lie  along  lines  of  curva¬ 
ture.  When  this  occurs,  boundary  conditions  c^n  be  approximated  at  best, 
and  then  only  with  the  introduction  of  extreme  mathematical  complexities. 
For  this  reason,  it  is  advantageous  lo  formulate  the  shell  theory  in  terms 
of  generalized  coordinates  so  that  boundaries  coincide  with  particulai 
coordinate  lines. 

This  section  presents  the  generalized  theory  upon  which  the  STAGS 
computer  code  is  based.  Although  no  attempt  is  madv.  to  be  exhaustive 
in  the  coverage  of  the  basic  shell  theory,  a  brief  description  of  the  funda¬ 
mental  aspects  is  given.  For  additional  material,  the  reader  is  referred 
to  Reference  9.  In  addition,  methods  for  computing  the  shell  middle  sur¬ 
face  input  parameters  are  presented. 

2.  1  Metric  of  the  Shell  Middle  Surface 

1  2 

Consider  a  surface  in  space,  described  by  coordinates  cp  and  cp  , 

which  is  embedded  in  a  three  dimensional  Euclidean  space  defined  by  the 

1  2  3 

Cartesian  coordinates,,  x  ,  x  ,  and  x  ,  as  shown  in  Figure  1.  The 
vector  r  to  any  point  on  the  surface  can  be  written  as 


where  the  are  unit  vectors  in  the  x  directions,  respectively.  Now 
consider  the  differential 


dr  =  dx  +  dx“  +  dxJ  k^ 


=  dx1  It. 


and  define  this  expression  for  dr  in  terms  of  the  shell  coordinates  cp 
and  cp2  b> 


dr  =  a^  dtp*'  +  a^  dcp2 


=  a  dtp 

a 


The  quantities  a^  ,  are  called  the  covariant  base  vectors  and  can  be 
written  as 

a  =  +  E  ^L.  +  k- 

a  1  hf  2  dcp"  3  d<p* 


1  Sep0' 


It  should  be  noted  that,  in  general,  the  base  vectors  a  are  not  unit 

a 

vectors  bu„  have  magnitudes  given  by 


Pi  •  al 


a2  ’  a2 


The  expression  for  incremental  arc  length  on  the  shell  middle  surface  is 


ds  =  dr  •  dr 
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(6) 


2  11  12  12 
ds  =  dcp  dcp  +  2  dtp  dtp  +  &2.Z  dtp  dtp 


affP  d/  d«pP 


The  quantities  a^g  are  called  the  components  of  the  covarimt  metric 
tensor,  and  are  defined  by 


a  a  =  a  •  a  „ 
ffP  ot  p 


Two  alternate  forms  of  Eq.  (6)  are 


?  211  12  222 
ds4  =  A  dcp  dcp  +2  C  dcp  dcp  +  B  dcp  dcp 


ds  ^  ^  dcp*  dcp*  +  2  A  B  cos  9  dcp*  dcp^  +  B^  dcp^  dtp** 


Both  of  these  formulations  have  been  used  in  the  STAGS  User's  Manuals; 
the  quantities  A,  B,  C,  and  0  are  related  to  the  components  of  the  metric 
tensor  by 


yf*ii 

Va22 


C  =  A  B  cos  9  -  a, 


cos  9  - 


all  a22 


It  can  be  seen  from  Eqs.,  '  7,  8)  that  A  and  B  are  measure0  m  the 

1  2 

length  along  the  coordinate  lines  cp  and  cp*-  ,  and  that  9  is  me  angle  between 


these  lines. 


Since  the  covariant  base  vectors  \  \  are  not  necessarily 

normal  to  one  another,  it  i*  sometimes  convenient  to  consider  a  set  of 
vectors  defined  by 


a 2  x  a3 


l 

-  i 


a3  x  al 


where  is  the  unit  normal  vector  as  shown  in  Figure  1  and  a  is  the 
determinant  of  the  metric  tensor 


al  X  a2 


all  a22  ‘  a12 


The  vectors  defined  by  Eqs.  (9)  ere  called  contravariant  base  vectors 
and  have  the  properties 


-1  _ 
a  •  a. 


-2  - 
a  •  a0 


a  •.  a„ 


-2  - 
a  •.  a. 


a  •  a 
o 


6 

6  is  the  Kr 
a 


•  kei  delta  and  nas  the  propertie1 


-s  - 


s 


fJ^TTV^VTS  rr?^4  T» *  ;<f^,^7<-avr>'rvr'^  CW' 


6e  .  (  1  '  * 
a  (0  ,  a 


The  contravariant  metric  tensor  components  are  defined  by 


oS  _  — o  —  P 
a  =  a  •  a 


which  can  be  written  in  terms  of  a  as 

a 


With  the  aid  of  Eqs.  (10  through  15),  Eqs,  (9)  can  be  written  in  the  form 


— O'  CV0  — 

a  =  a  ar 


2.  Z  Curvature  Tensor  of  the  Shell  Middle  Surface 


The  curvature  and  twist  of  a  surface  are  defined  by  the  curvature 
and  torsion  of  lines  embedded  in  the  surface  relative  to  the  unn  -ormal  vec¬ 
tor.  For  instance,  the  normal  curvatures  of  a  surface  with  respect  to  the 
1  2 

coordinate  lines  cp  and  cp  are  defined  by 


bl  '  a3 


b2  '  a3 


whore  a.  -u.d  b  re  the  .  •  ’’  t  .ri  vo.  T,,r&  .it  t':e  >  ooru.  •  s  1 1  u'.rs 


and  cp  ,  respectively.  These  vectors  can  be  written  as 


where  ds^  and  ds^  are  given  by 

dsl  =  v/all  d:pl 
ds2  =  ^11  d^ 


Hence,  the  normal  curvatures  of  the  shell  middle  surface  are 


1  d 

dcp' 


1 

'/a22 


(17) 


The  twist  of  a  surface  with  respect  to  a  coordinate  line  is  the  tor¬ 
sion  of  the  coordinate  line  wi>.h  the  sign  chosen  such  that  a  positive  twist 
occurs  when  the  normal  vector  a^  rotates  about  one  coordinate  line  towards 
the  other  coordinate  line.  This  leads  to  the  definitions  for  twist  of  a  surface 


b 


tl 


b 


t2 


'/SU 


(18) 
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The  curvature  tensor  of  a  surface  is  defined  by 


(19) 


This  is  a  symmetric  tenso’*  whose  indices  are  raised  and  lowered  accord¬ 
ing  to 


a’* 

apP 


(20) 


The  components  of  b^g 
and  twist  by 


are  related  to  the  changes  iu  normal  curvature 


b 

no 


aat 


/a— 

aa 


(o  not  summed) 


■(at  P;  a,  P  not  summed) 


(21) 


It  can  be  seen  with  the  aid  of  £qs.  (10,  15,  20)  that,  for  orthogonal  coordi¬ 
nates  (a12  =  0),  b  ^  =  bt2  , 

Two  invariants  associated  with  the  curvature  are  the  mean  curva¬ 
ture  H  and  the  Gaussian  curvature  K 


H  s 


1 

I 


ha 

a 


K 


b 


a 

p 


) 


(22) 
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2.  3  '  Displacements  and  Derivatives 

The  displacement  vector  of  a  point  on  the  shell  middle  surface  is 
defined  by 


u  =  a  +  w  a_ 
or  3 


(23a) 


ry 

The  quantities  u  are  called  the  contravariant  components  of  in-plane 
displacement.  An  alternate  form  of  (23a)  is 


—  a  — 

u  a  +  w  a, 
o  3 


(23b) 


where  the  u  are  the  covariant  components  of  in-plane  displacement. 

Since  a,  is  a  unit  vector  and  normal  to  a  and  a®  ,  there  is  no  dis- 
3  0- 

tinction  between  covariance  and  contravariance  for  w  .  From  Eq.  (16) 

CK 

it  can  be  seen  that  u  and  u  are  related  by 

a 


a  _  oP 
u  =  a  u 


u  =  a  R  u 
O'  OP 


(24) 


The  covariant  and  contravariant  displacement  components  can  be 
related  to  physical  quantities.  Consider  the  displa  ce'  ient  vector  written 
in  terms  of  unit  covariant  vectors 


^Ti 


/r.; 


22 


w 


(25) 


The  quantities  u  and  v  represent  the  physical  components  of  in-plane 
displacement  xn  the  directions  defined  by  the  .covariant  base  vectors  aj, 
and  a^  ,  respectively.  Equating  Eqs.  (22.  to  Eq.  (25)  with  help  from 
Eqs.  (7,  14,  16)  yields 


t  l\a  l\",  im-  i  l  '  iJntVit  i 


°l  =  “ 


_  *12  , - 

u-  -  -  u  +  /a00  v 

2 


The  partial  derivative  of  the  covariant  base  vectors  are  defined 


Fq  a  +  b  _  a, 
Per  P  a  B  3 


The  b  -  a^  term  can  be  deduced  from  Er.  (19).  It  follows  from  Eq.  (In' 


that  T  o  can  be  written  as 
ffP 


rP  =  rp  -  a  tP 


17  ■  a 


The  partial  derivatives  of  the  contravariant  base  vectors  can  be  deduced 
from  Eq.  (28,  29)  as 


r~a  —  p  ,  a  — 

a  f  o  -  a  0 
£p  fc  3 


The  quantities  7  ~  are  Christoffel  symbols  of  lb 


?cond  kind.  Symbols 


1 


of  the  first  kind  are  defined  by  lowering  the  upper  index 


*arPp  ^oP  a 


*XP  1^"  ap 


These  quantities  expressed  in  terms  of  the  partial  derivatives  of  a 


H 


3aD  9a  a\ 

Pp  _  crp 


With  the  aid  of  Eqs.  {28,  30),  tbe  partial  derivative  of  the  displace¬ 
ment  vector  u  [Eq.  (22,23)]  can  be  written  as 


— -  =  uP|  aft  +  uf  a^  -  b^waft  +  a 

cuf  P  O'  P  3  or  P  aq)or  3 


I  —  P  ,  -  P  -  ,  p  —  dw  — 

=  u  aK  +  b  uaa,  -  b  w  aQ  +  -  a~ 

P'o  »  P  3  O'  P  ^  o'  3 


The  quantities  Upj^  and  up|^  are  called  the  covariant  derivatives  of 
Up  and  uP  with  respect  to  ya  and  are  given  by 


0.1  -  -i  -  r?  u 

p  »  ZJX  Pa  P 


upi  »  +  rs  uO 

“  9/  p« 


The  concept  of  covariant  differentiation  can  be  extended  to  second  order 
tensors  such  as  b  _  . 


db 


3cp 


£§  -  rx  b 

op  0X 


bcA 


a,* 

-°L  + 

dcpP 


r\  bp 
pX  o 


(35) 


The  quantities  appearing  in  Eqs.  (33)  can  be  regrouped  to  define  the  dis¬ 
placement  gradients  and  p^ 


M  ■  Ve  ‘  b«P  w 


„  3w  .  ,  P 

P  =  -  +  b H  u 

"  Scp°  "  P 


(36) 


Hence,  the  derivatives  of  u  take  the  form 


a 


(37) 


2.  4  Deformation 

The  deformation  of  the  shell  middle  surface  can  be  specified  in 
terms  of  the  changes  in  the  metric  and  curvature  tensors.  With  the  de¬ 
formed  state  of  the  shell  characterized  by  a  tilda  (a  .  ,  £T  )  ,  the 

op  op 

strain  and  curvature-change  tensors  can  be  defined  by 


2  £oP  "  aaP  "  a<*0 


a 


a  a 


(38) 
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It  can  be  seen  with  the  aid  of  Eqs.  (6,  13,  22)  that  these  definitions  lead 
directly  to  the  changes  in  incremental  arc  length  and  mean  curvature 


ds^  -  ds^  =  2  e  .  dtp0,  dcp^ 


H  -  H  = 


1  K-  B  » « 

2  o’  p 


Although  e  and  k  p  are  independent  of  any  particular  metric  tensor, 
op  o 

it  is  convenient  to  refer  these  quantities  to  the  metric  of  the  undeformed 

middle  surface;'  i.  e.  ,  in  operations  expressing,  for  instance,  covariant 

strain  or  curvature-change  tensors  in  terms  of  contravariant  tenrors  (or 

vice  versa)  a°^  and  a  .  are  used  rather  than  aa ^  and  a  .  „  For 

op  op 

example, 


e*  =  c 
o  op 


•  P 

h  -  n  a  a 
OP  O  p  p 


In  addition  to  strain  and  curvatur e-change,  portions  of  the  shell  middle 

surface  may  undergo  finite  relations.  If  such  is  the  case,  the  expressions 

for  e  n  and  h  „  ,  when  written  in  terms  of  the  displacement  gradients, 

OP  OP  re. 

must  reflect  this.  Since  the  general  expressions  for  e  „  and  n  _  are 

b  r  op  a  P 

extremely  complicated,  it  is  desirable  to  use  simpler,  approximate  expres¬ 
sions  whenever  possible. 

The  rotation  of  any  part  of  the  shell  middle  surface  can  be  split 


'  n  is  not  generally  symmetric;  therefore,  its  indices  must  be  moved  straight 
up  and  down  only.  The  dot  appearing  over  the  o  assures  this..  Note  that  a 


dot  m  the  terms 


is  not  required  since  e 


e  ,  ;  hence,  for 

Pa 


symmetric  tensors,  no  distinction  m  the  ordering  ol  indices  is  mjmri  d. 
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into  two  parts,  an  out-of -plane  rotation  or  "tilt"  and  an  in-plane  rotation 
commonly  called  the  "rotation  about  the  normal.  "  When  the  angle  of  tilt 
(fl)  is  moderate,  the  tilt  and  in -plane  rotation  (u>)  can  be  approximated  by 


sin  <c  a- 


Tjz  (Y21 ' Vl2) 


(41) 


sin  ^  Cl  «  pQ  P 


a 


The  expressions  for  e  „  and  h  .  used  in  STAGS  are 

op  Op 


=  I  (Ya'fi  +  Yfta)  +  2  YP*  ''■’’p  +  I  S»  PP 


Ka  =  Pr  !  +  bP  y  q  -  b  v P 

o P  Ptt  O'  TP P  Op  p 


(42) 


These  approximations  are  based  upon  the  assumptions  that  the  tilt  earn  be 
moderately  large  (Q  <  .  3)  and  that  the  in-plane  rotation  is  of  the  same 
order  of  magnitude  as  the  square  root  of  a  typical  middle  surface  strain 
(ui  =  0  (/e)  }.  A  complete  derivation  of  the  above  is  given  in  Reference  9. 
Physical  components  of  strain  and  curvature-change  for  lines  of  curvature 
coordinates  are  given  in  Section  2.  9. 

2.  5  Strain  Energy 

The  strain  energy  density  for  thin  elastic  shells  is 


U 


1 

1 


E 

- -j 

« 

1- V 


(1-v) 


ap  PX 


4-  v  a 


oP  p  X 


(43) 


For  shell  coordinates  §  and  T)  and  with  the  use  of  Eqs.  (8,  15),  Eq.  (43) 


can  be  cast  in  the  form 


.  -.-4  2  4A  cos  9 


n  =  2“  j  (A  Sin  9)  -  - g -  (A  sin  9) 


+  2  [l-(l-v)  sin29]  (AB  sin29)'2  e?§  e  ^  +  2  (AB  sin29)-2  [(1-v)  +  (1+v)  cos2 9]  e2^ 


^TT1  iB  si"  6>'4  «5n  eT|T|  +  <E  8in  e>'4  4 


K  )  /A  .  „«-4  2  4A  cos  9  /A  .  -4 

+  T  )  (A  Sln  0)  k5§  "  — b -  (A  sin  8/  H§?;  h|ti 


+  2  [l-(l-v)  sin2  8)  ]~2  (AB  sin29)~2  +  2  (AB  sin2  8)  2  [(1-v)  cos29] 


4B  cos  9  lD  .  „.-4  ,  /rj  ..-4  2  I 

- 5—  (B  sm  6!  +  (B  sin  6) 


where  D  and  K  are  the  membrane  and  bending  stiffnesses,  respectively. 


Et 
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The  covariant  components  of  strain  and  curvature-change  expressed 
in  terms  of  the  displacement  gradients  are 


I 

I 
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*55  Ys§  +  r  (Y5§Y-S  + 

*SH  =  I  (Y5H  +  YT)5  +  Y§§  V*%  +  yT\ly\  *  (46) 

®T\T)  =  YT|T1  +  2  (Y?TlY-\  +  YT1T1  y.\  +  ?T1  N) 
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where  the  commas  denote  partial  differentiation.  The  displacement  gradients 
written  as  functions  of  the  physical  components  of  displacement  u  and  v  are 
(see  Eqs.  (8,  26,  27,  34,  36)) 


B5  3  w,?  +  K  b55“  *  W  b 5T1  v 


Bt)  =  W,H  +  X  b5Ti  u  +  W  blffl  v 


(48) 
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Ac.,  ^  +  A/B  {A,  ^  -  B,  ^  cos  0)  v  +  A  cos  0  v,  ^  w 
Au,  ^  +  [(A  cos  0),  ^  -  B,  g  ]  v  + Aios  6  v,  -  b^w 
[(B  cos  8),  |  -  A,  ^  ]  u  +  B  cos  6  u,  ^  +  B  v,  ^  ■:  b^w 
B/A  (B,^  -  A,  cos  0)  u  +  B  cos  B  u,  ^  +  B  v,  ^  w 

1/A  u,  ^  +  [(cos  0/  (AB  sin2  0)  ]  [A,  ~  (B  cos  9),  ^  ]  u 


+  [1/  (AB  sin2  0)  ]  [A,  ^  -  cos  0  B,  ^  ]  v  -  b^  w 

1/A  u,  ^  +  [1/ (A2  sin20)][(A  cos  0),  ^  -  B,  ^  ]  v 

2  2  E 

+  [cos  0/  (A  sin  0)]  [A,  ^  cos  0  -  B,  F  ]  u  -  b!^  w 

1/B  v,  ^  +  [1/(B2  sin2  6)  j  [(B  cos  0),  -  -  A,  ]  u 

+  [cos  0/ (B2  sin20)  ]  [B,  ^  cos  0  -  A,  ^  ]  v  -  b^  w 

1/B  v,  +  [cos  9/ (AB  sin2  6)1  [B,  *  -  (A  cos  9),  _]  v 

+  [i/(AB  sm20)]  [B.  ?  -  cos  0  A,  ^  1  u  -  b^1  w 

■S  '  U 


(49) 
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and  the  Christoffel  symbols  are  (Eqs.  (8,  32)) 


4  ■ 

[BA,  ^  +  (AA,  ^  -  C,  cos  8]/(AB  sin2  8) 

‘§5  " 

[C,  ^  -  BA,  ^  cos  0  -  AA,  ]  /(B2  sin2  8) 

r§  = 

1  5T) 

[A,  ^  -  B,  ^  cos  0]  / (A  sin2 8) 

rl  - 

[B,  -  -  A,  cos  0]  / (B  sin2 8) 

5  h 

r§  = 
Tin 

[C,  ^  -  BE,  ^  -  AB,  cos  0]  /(A2  sin20) 

r  ^  = 

T|I| 

[AB,  _+  (BB,  .  -  C,  ~)  cos  0 J  /  (AB  sin20) 

Substitution  of  Eqs.  (48-50)  into  Eqs.  (46,47)  yields  the  covariant 
strain  and  curvature-change  tensors  as  functions  of  the  physical  displace¬ 
ment  components  u  ,  v  ,  and  w.  These  equations  are  then  substituted 
into  Eq.  (44)  to  obtain  the  strain  energy  as  a  function  of  the  displacements. 


The  effects  of  geometric  imperfections  have  been  accounted  for 
by  modifying  Eqs.  (46)  to  include  small  values  of  an  initial  normal  dis¬ 
placement  w  .  The  terms  w,  _  (3,  ,  w,  _  0..  and  w,  _  0_  +  w,  „  0_  were 

SIMM  l  5  h  5 


added  to  the  three  middle  surface  strain  eeir  ,  e„_  ,  and  ecr, 

5s  Tim  s'I 


respectively. 

Geometric  imperfections  are  important  because  the  critical  loads 


of  many  shells  are  sensitive  to  such  imperfections.  In  addition,  there  ari 
many  cases  where  there  exist  planes  of  symmetry  with  respect  to  loading 
and  geometry.  In  such  cases,  antimetric  deformations  will  only  be  found 
if  they  are  "triggered"  by  the  inclusion  of  antimetric  geometric  imper¬ 
fections. 


2.  6  Solution  Procedure 

The  solution  procedure  used  in  STAGS  is  based  on  the  principal  of 
stationary  potential  energy.  After  the  expression  for  strain  energy  den¬ 
sity  is  formed,  as  explained  in  the  previous  section,  the  displacements 
and  their  derivatives  are  replaced  by  appropriate  finite  difference 
expressions.  {A  set  of  finite  difference  expressions  for  grids  with 
variable  spacing  is  described  in  Section  3.  0. )  The  strain  energy  density 
at  mesh  station  i  is  then  written  in  the  form 


AUl  =  j  Zl*  D1  Z1 


where  D  is  a  6  x  6  matrix  of  constants  and  Z1  is  the  vector  of  strains  and 
curvature  changes  at  station  i  .  (In  this  report  all  vectors  are  understood 

i# 

to  be  column  vectors  and  *  desig.iates  the  adjoint  operator.  Thus,  Z 
is  a  row  vector. )  The  matrix  D*  is  obtained  by  integration  through  the 
shell  wall  and  is  a  function  of  the  geometric  parameters  of  the  shell  and 
of  the  material  properties.  The  strain  vector  Z1  is  a  nonlinear  (quadratic) 
function  of  the  displacement  unknowns  and  the  geometric  parameters.  The 
vector  of  stress  resultants  at  station  i  is  given  by 


Sl  =  D1  Z1 


The  total  strain  energy  of  the  shell  is  then 


U  = 


E  4Ui 


i  th 

where  a  is  the  area  of  the  i  mesh  subregion.  The  potential  energy  of 
the.  work  done  by  external  forces,  W  ,  may  be  expressed  in  discrete  form 
by 


W  =  X';:  F 
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where  F  is  the  vector  of  external  forces.  As  the  strain  expressions 
are  of  second  order  in  the  displacement  components,  the  total  potential 
energy,  V  ,  of  the  shell  is  a  polynomial  of  4^  degree  in  the  discrete 
displacement  unknowns.  It  is  given  by 

V  =  U  -  W  (55) 

A  necessary  condition  for  static  equilibrium  is  that  the  potential  energy  be 
stationary.  For  a  polynomial,  this  requires  that  the  gradient  of  V  vanishes 
and  leads  to  the  equation 

LX  =  F  (56) 


Here  L  is  defined  as  the  nonlinear  operator  such  that 

LX  =  Grad  U  (57) 

The  derivation  of  the  complete  nonlinear  solution  of  Eq.  (56)  as  well  as  of 
bifurcation  buckling  is  facilitated  by  introduction  of  the  concept  of  the 
derivative  L'  of  L  (Ref.  10).  In  particular,  for  the  operator  L  ,  the  deri¬ 
vative  L'  (sometimes  called  the  Frechet  derivative  of  L)  is  an  n-by-n 
matrix  whose  elements  are 


S2U 


axii>  axo) 


(58) 


Because  L'  is  a  function  of  a  particular  displacement  vector  X  (unless 
the  nonlinear  terms  are  dropped),  the  Frechet  derivative  will  usually  be 
denoted  L'  to  indicate  this  dependence.  With  the  use  of  the  derivative 
L'  of  the  operator  L  ,  Newton's  method  may  be  readily  generalized  to 
obtain  the  solution  of  Eq.  (56),  The.  iteration  is  defined  by 

Sf  <xk+i  -  xk>  -  F  -  LXk 


(59) 


If  Xq  is  an  estimate  sufficiently  close  to  the  solution  X  and  if  Lj^  is 
not  a  singular  matrix,  die  iteration  converges  to  X  .  Under  these  assump¬ 
tions,  it  also  can  be  shown  that  the  converged  solution  is  unique  (Ref.  10). 

2.  7  Bifurcation  Buckling 

The  application  of  Newton's  method  and  the  modified  Newton  method 
in  STAGS  to  obtain  a  nonlinear  collapse  analysis  is  discussed  in  the  pre¬ 
vious  section.  Jt  is  interesting  to  note  that  the  mathematical  characteriza¬ 
tion  of  bifurcation  buckling  also  is  provided  by  the  generalized  Newton 
method.  Let  XQ  be  a  solution  of  Eq.  (56)  under  a  given  vector  F  of 
external  forces.  If  every  neighborhood  of  Xq  contains  another  vector  Y 
which  satisfies  the  equation 

LY  =  F  (60) 

then  bifurcation  is  said  to  occur  for  the  shell  under  the  load  F  .  From  the 
previous  remarks  on  the  conditions  for  convergence  of  Newton's  method  to 
a  unique  solution,  it  follows  that  a  necessary  condition  for  bifurcation  is 
that  ^Xq  be  a  singular  matrix,  or 

det  (L^0)  =  0  (61) 


Classical  bifurcation  buckling  theory  (with  linear  prebuckling  analysis.)  may 
be  easily  obtained  from  Eq.  (61).  It  is  assumed  that  Xq  may  be  written 

XQ  -  XXL  (62) 

where  is  the  linear  solution  for  a  lo?  d  vector  F^  .  Thus,  Eq.  (61) 

becomes 


det 


=  0 


(63) 
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Equation  (63)  is  an  algebraic  eigenvalue  problem  of  the  form 


det  (A  -  XB  -  X2C)  =  0  (64) 

In  classical  bifurcation  theory,  the  C  matrix,  which  arises  from  the  pre¬ 
buckling  rotations,  is  often  omitted  and  the  eigenvalue  problem 

AX  =  XBX  (65) 


is  obtained. 


When  bifurcation  is  considered  but  the  prebuckling  displacements 
are  not  linear,  the  solution  of  Eq.  (61)  generally  requires  a  stepwise  pro¬ 
cedure.  One  such  method  is  given  by  the  recurrence  equations 


det 


=  0 


X,  =  h  .  X, 

krl  k+1  k 


(66) 


in  which  the  starting  vector  Xq  may  be  represented  by  the  linear  solution. 

A  sequence  of  eigenvalue  problems  is  solved  and,  if  the  method  is 
successful,  ^  approaches  unity.  A  nonlinear  bifurcation  treatment  [equiva¬ 
lent  to  Eq.  (66)]  was  presented  in  Reference  11  and  has  been  used  successfully 
to  study  a  large  variety  of  problems.  Ecr  the  two-dimensional  problems 
under  consideration  here,  it  appears  that  such  methods  may  be  as  costly  as 
the  complete  nonlinear  analysis  available  in  STAGS.  Consequently,  only  a 
classical  bifurcation  buckling  analysis  is  implemented  in  the  STAGS  pro¬ 
gram. 

The  formation  of  the  A  and  B  matrices  of  Eq.  (65)  will  be  con¬ 
sidered  briefly.  The  e'ements  of  the  Frechet  derivative  matrix 
(which  define  the  natrices  A  and  B)  are  determined  according  to  Eq.  (58). 
The  rules  for  computing  derivatives  of  polynomials  are  easily  programmed, 


and  the  formation  of  the  a  and  B  matrices,  therefore,  is  we?l  suited  to 
automatic  treatment  on  the  computer.  Thux,  for  example,  if  and 

X^j,  are  the  i^1  and  jfb  displacement  components,  we  have,  using  Eqs. 
(52),  (53),  and  (54): 
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Examining  the  k^  term  of  this  sum, 


a2z.uk 

^1X^ 


s2zk?; 


2X(i)  2X(j) 
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(67) 


(68) 


In  the  first  term  on  the  right-hand  side  of  Eq.  (68),  note  that  S  is  the 

linear  stress  resultant  vector  at  station  k  and  that  only  the  quadratic 

terms  (rotations)  need  be  considered  in  forming  the  partial  derivatives 
2  k* 

a  Z  /3X  3X  .  Contribute  s  from  this  term  go  into  the  B  matrix. 
Assuming  the  prebuckling  rotations  may  be  neglected  for  the  classical 
theory,  the  last  term  of  Eq.  (68)  generates  additions  only  to  the  A  matrix. 
The  A  matrix  is  then  identical  to  the  linear  stiffness  matrix.  If  the  pre¬ 
buckling  rotations  are  included  (nonlinear  bifurcation),  the  last  term  of 
Eq.  (68)  generates  a  C  matrix  and  provides  additional  contributions  to  the 
B  matrix.  In  this  case,  the  prebuckling  stress  resultant  vector  S  would 
also  include  nonlinear  terms. 

In  conclusion,  it  should  be  noted  that  bifurcation  buckling  theory  is 
often  based  on  the  concept  of  adjacent  equilibrium  states.  Of  course,  the 
same  algebraic  eigenvalue  problem  is  ultimately  obtained  by  both  methods. 
However,  the  approach  presented  here  seems  to  provide  a  more  simple 
recipe  for  definition  of  the  basic  matrices  of  Eq.  (65).  The  recipe  is  out¬ 
lined  m  Eqs.  (67)  and  (68)  and  leads  to  str aightforward  algebraic  procedures. 
In  addition,  the  relations  between  linear  and  nonlinear  bifurcation  theory 
and  Newton's  method  are.  clarified. 


2.  8  Computation  of  the  Components  of  the  Metric  and  Curvature  Tensors 

The  components  of  the  metric  tensor  A,  B,  C  and  of  the  curvature 

tensor  b  „  can  be  computed  according  to  the  following  sequence  of  opera- 
Of  P 

tions: 

1)  Denote  the  cartesian  coordinates  x*  ,  x2  ,  x2  by  x,  y,  z  , 

respectively,  and  the  surface  coordinates  ©  ,  cp2  by  §  ,  Tj .  The  x  ,  y  , 
z  are  then  determined  as  functions  of  §  and  T|  . 

x  =  x  (§  ,  Tj) 

y  =  y  (I ,  T)  (69) 

z  =  2  (5  ,  Tj) 


2)  The  metric  tensor  components  are  computed  with  the  aid  of  Eqs.  (4,  7) 
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The  components  of  the  curvature  tensor  are 
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where 
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Since  the  sign  of  the  curvature  tensor  was  chosen  so  that  positive 
curvature  results  from  an  inward  point  unit  normal  [see  Eq.  (19)3  >  care 
must  be  taken  to  preserve  the  sense  of  b  ^  by  remembering  that  it  is 
defined  by  the  cross  product  a^  x  a^  . 

As  an  example,  consider  the  case  of  an  elliptic  cone  as  shown  in 
Figure  2.  The  parameters  o  and  (3  are  the  tangents  of  the  cone  half  apex 
angles  in  the  x-z  and  y-z  planes,  respectively,  §  is  the  elliptic  coordinate, 
and  7]  is  the  axial  coordinate.  The  relationships  between  x,  y,  z  and 
l  ,  T1  are 

x  =  o T|  cos  § 

y  =  PT|  sin  §  (73) 

z  =  T1 

Note  that  this  choice  of  §  and  T]  results  in  an  outward  pointing  normal  as 
shown  in  Figure  2. 

The  a  D  and  b  are  computed  from  Eqs.  (70,  71) 
op  op 
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2.  9  Physical  Components  of  Strain  and  Curvature-Change 

Fo-  lines  of  curvature  coordinates,  the  physical  components  of 
strain  and  curvature-change  are  given  by 

Vei  =  7===-  <nosum) 

1  aa 


(no  sum) 
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With  the  components  of  the  metric  and  curvature  tensors  written  as 
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The  physical  components  of  the  displacement  gradients  (Eq.  (36))  take  the 
form 
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V)  =  w,i  +  57  “a) 
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where  a  comma  denotes  partial  differentiation  and  the  quantities  are 
the  physical  components  of  displacement  in  the  directions. 

The  physical  components  of  strain  and  curvature-change  are  deter¬ 
mined  from  Eqs.  (42)  and  are  found  to  be 


(11) 

Y(H)  +  2 

v(ii)  + 

1  j 

(12) 

2  ( Y(12)  4 

Y(21)) 

1 

2 

(22) 

■(22;  +  2 

Y(12) 

)  +  i 


,  1  q2 

Y  T 


Y(ll)  \l2)  +  2  Y(2i;  y(22)  +  2  p(i)  "(2) 


(78) 


I  2 
T  V, 


(22) 


?(2) 


Si 


Section  3.  0 

FINITE  DIFFERENCE  GRID  WITH  VARIABLE  SPACING 


For  better  economy  in  the  analysis  a  capability  has  been  provided 
for  the  use  of  variaole  .spacing  finite  difference  grids.  The  shell  surface 
is  covered  with  a  system  of  mesh  lines  (see  Fig.  3)  whose  coordinates 
are  given  by 


x.  ,  i  =  1,  m 

l 

and 


Yj  »  j  =  1«  n 

where  x  and  y  are  the  curvilinear  surface  coordinates.  Corresponding 

to  each  pair  of  values  (i,  j)  i  =  1,  m;  j  =  1,  n,  a  rectangular  region  R.  . 

L  J 

is  defined  with  sides  of  length 


a.  , 
i.  J 


1/2  h+i  -  xi-i  i 


b.j  =  1/2  'Vi  '  yj-i 1 

Note  that  R^  ^  is  rectangular  on  the  map  of  the  shell  provided  by  the  surface 
coordinates  but  not  generally  on  the  shell  itself. 

The  regions  R.  (and  lengths  a.  .  ,  b.  .)  are  modified  at  bounda- 
L  J  J  J 

ries  of  a  shell  by  including  only  those  portions  which  are  msicie  the  panel. 

A  double  integral  of  a  function  f  over  the  region  R  of  the  panel  may  then 
be  approximated  by 
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(80) 
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The  discretization  is  completed  when  the  integrand  functions  f.  ^ 

are  evaluated  at  the  centroids  o:  the  regions  R.  .  in  terms  of  the  neighbor- 

"  i*  J 

ing  displacement  components.  It  shou^  first  be  noted  that  the  tangential 

displacements  u  and  v  have  been  located  at  corners  of  the  regions  R.  .  . 

J 

Furthermore,  the  energy  expressions  for  a  shell  include  derivatives  of  u 
and  v  only  up  to  the  first  order.  Hence,  even  with  arbitrary  spacing,  only 
central  difference  formulas  for  the  u  and  v  displacements  are  required. 

In  contrast,  the  normal  displacement  w  has  been  located  at  the  mesh  node 
points  (x^  ,  y^)  and  more  general  finite  difference  formulas  must  be  developed. 

The  coordinates  of  the  centroid  '/  a  region  R.  .  are  given  by 

J 


x.  =  1/4  (X..J  +  2x.  +  *.+1) 

y.  =  1/4  <yj.I  +  Zy;  +  yjtl) 


(■31) 


Variable  spacing  is  first  considered  with  respect  to  a  single  coordinate  x 
only.  With  the  help  of  a  Taylor's  expansion  (or  equivalently  by  the  use  of 
interpolation  formulas),  the  difference  formulas  for  w  ,  w,  x  and  w, 
at  may  be  established  as 

(w).  =  w|~  =  w._1/16  •  [(h-k)  •  (3k+h)/(h2+  hk)] 

i 

+  w7l6  [(h+3k)  •  (3h+k)/(h+k)]  (82) 

+  w  +1/ 16  •  [(k-h)  •  (Jh+k)/(hk+k2)] 


(w,x)  s  w,xi_  =  -  w  /(2h) 

1  v  ‘ 


X. 

1 


+  w.  [1/  (2h)  -  1/  (2k)  ] 


wi+i/ (2k) 


(83) 


xx\  ~  w’  XX 


X. 

1 


=  ^.1  •  2/  [h  >  (h+k)] 


•  2/  (h  •  k) 


(84) 


+  w.+1-  2/  [k  •  (h+k)] 


where 


h  =  x.  -  x.  . 
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The  corresponding  formulas  for  the  y  coordinate  are  obtained  by 
appropriate  substitutions  and  are  denoted  with  superscripts 
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The  required  two  dimensional  diTeroice  formulas  are:  now  obtained 
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by  combining  the  formulas  for  both  coordinate  directions 


(87) 


In  general,  these  equations  involve  the  9  point  "star"  of  neighboring  values. 
However,  it  is  easilv  seen  that  all  of  the  formulas  reduce  to  the  standard 
central  difference  formulas  when  uniform  rectangular  spacing  is  used.  All 
of  the  difference  formulas  are  exact  when  the  displacement  function  w 
behaves  quadratically 

The  inclusion  of  nonorthogonal  coordinates  and  of  variable  grid 
spacing  extends  considerably  the  class  of  cases  that  can  be  solved  by 
use  of  STAGS.  The  grid  lines  can  be  made  to  follow  boundary  lines  and 
cutout  edges  rather  than  lines  of  curvature  on  the  shell  surface.  As  an 
example,  for  a  cylinder  with  a  circular  cutout,  one  can  use  a  grid  as 
shown  in  Figure  4.  This  grid  system  is  described  by  the  mapping 
function 

f(x)  =  a(l-x)  + 

g  (x)  -  b(l-xj  +  Lx 
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Grid  for  Cylinder  With  Circular  Cutont 
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X  =  R  sin  _F  cos  0/ Rl 


(88) 


Y  =  R  cos  If  cos  0/R  ^ 

Z  =  F  sin  0 

However,  it  was  found  that  use  of  this  grid  leads  to  inaccurate  results 
unless  the  spacing  between  the  straight  gridlines  are  very  close  together 
in  the  neighborhood  of  the  corner.  The  reason  for  this  is  that  the  equations 
for  the  strains  are  inaccurate  if  the  angle  bet-.een  the  coordinates  changes 
significantly  between  two  adjacent  gridpoints.  It  appears  to  be  more 
practical,  until  other  finite  difference  expressions  have  been  developed, 
to  use  a  different  approach.  For  shells  with  cutouts  that  cannot  easily  be 
made  to  follow  a  regular  net,  it  is  suggested  that  the  user  written  subrou¬ 
tine  for  a  variable  thickness  shell  be  used.  The  shell  thickness  in  set 
equal  to  zero  if  a  gridpoint  falls  inside  the  cutout.  The  computer  program 
then  eliminates  as  unknowns  th'*  displacement  components  at  points  with 
zero  thickness.  This  method  is  demonstrated  in  the  example  given  in 
Section  5.  3. 


Section  4.  0 

INELASTIC  BEHAVIOR 


4. 1  Introduction 

Due  to  the  extreme  complexity  of  the  problem,  it  has  been  necessary 
to  fcrmulate  theories  of  plasticity  which  greatly  simplify  material  behavior. 
While  in  many  cases  these  theories  give  satisfactory  results,  theie  are 
other  cases  in  which  they  fail.  It  is  shown,  for  instance,  in  Reference  12 
that  for  loading  histories  with  sharp  turns  in  the  stress  space  the  classi¬ 
cal  theory  with  isotropic  strain  hardening  may  give  very  poor  results. 
Typically  at  collapse  there  is  a  very  sharp  change  in  deformation  pattern 
and,  consequently,  a  sharp  turn  in  the  stress  path.  Other  theories  have 
been  proposed  which  more  adequately  describe  the  material  behavior  in 
such  cases  than  does  the  classical  theory.  The  Batdorf-Budiansky  slip 
theory  (Ref.  13)  is  probably  too  cumbersome  for  practical  application,  but 
the  type  of  theory  proposed  by  White  (Ref.  14)  and  Besseling  (Ref.  15) 
appears  very  promising  because  it  is  rather  simple  in  its  application,  yet 
it  retains  such  features  as  strain  hardening  and  the  Bauschinger  effect. 

For  these  reasons,  it  was  selected  for  use  in  the  STAGS  code. 

Introduction  of  inelastic  behavior  has  been  done  within  the  frame¬ 
work  of  the  energy  principle  upon  which  the  elastic  analysis  was  based. 
Essentially,  the  plastic  deformations  are  considered  as  load  terms;  they 
are  completely  analogous  to  thermal  expansions  except  that  they  are  not 
known  in  advance.  A  series  of  clastic  problems  are  solved  by  the  use  of 
energy  principles  in  which  the  "load  terms"  are  gradually  modified  until 
they  correspond  to  "he  computed  state  of  stress  and  to  specified  nonlinear 
stress  strain  relations  at  all  points  over  the  shell  surface,  and  through  the 
shell  thickness. 


4.  2 


The  White -Besseling  Theory 


The  tfhite-Besseling  Theory,  as  applied  here,  assumes  that  the 
material  consists  of  r  everal  components  which  have  identical  elastic  pro¬ 
perties  and  exhibit  ideal  plasticity  (no  strain  hardening)  but  have  different 
yield  strength.  As  the  strain  is  the  same  in  all  components,  the  stress- 
strain  curve  will  experience  a  decrease  in  slope  as  the  stress  reaches 
the  yield  limit  for  any  one  of  the  components;  the  respective  components 
then  cease  to  take  additional  load.  The  composite  thus  exhibits  strain 
hardening  with  a  piecewise  linear  stress -strain  relation.  Use  of  only  one 
component  will,  of  course,  result  in  application  of  ideal  plasticity  theory. 

If  the  stress  is  reversed  after  loading  beyond  the  yield  limit  for  one  or  more 
components,  yield  will  occur  in  the  reversed  direction  at  an  average  stress 
in  the  composite  which  is  lower  than  the  stress  for  original  yield.  This 
is  demonstrated  in  the  uniaxial  stress-strain  curve  shown  in  Figure  5. 
Tension  is  first  applied,  OAB,  beyond  the  yield  limit  and  followed  by 
strain  reversal,  BCD,  into  the  zone  of  yield  in  compression.  The  yield 
ellipse  for  tne  weakest  component  and  the  loading  history  in  this  component 
are  also  shown  in  this  figure.  Clearly,  yield  in  compression  will  occur 
when  the  total  strain  is  (e^  ~  2e  ),  i.  e.  ,  the  yield  in  compression  occurs 
at  a  much  lower  stress  if  the  material  previously  has  been  subjected  to 
tension  stresses  above  the  yield  point.  To  introduce  the  Bauschinger 
effect  this  way  is  appealing  because  it  reflects  the.  microstress  theory 
which  now  generally  is  accepted  as  the  explanation  of  the  Bauschinger  effect. 

4.  3  Implementation  of  the  White-B essehng  Theory  in  STAGS 

The  White-Bessehng  plasticity  theory  is  implemented  in  the  com¬ 
puter  program  in  the  following  manner: 

1)  The  inelastic  behavior  of  the  material  is  defined  through 
specification  of  the  uniaxial  stress-strain  curve.  This  curve  is  piecewise 
linear  and  the  input  consists  of  stress  and  strain  at  each  of  its  comers.  The 
relative  volume  anc  the  yield  strength  for  e..ch  of  the  components  is  then 
computed  internal!  .  The  transverse  strain  is  determined  such  that  the 
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Figure  5  Stress  Strain  Curve  Showing  Bauschinger  Effect 


total  stress  in  the  transverse  direction  is  zero.  It  seems  possible  to 
reflect  more  accurately  the  actual  stress -strain  relations,  including 
anisotropy,  if  self-equilibrating  initial  stresses  are  included. 

(2)  The  strains  are  estimated  for  all  points  in  the  shell  over  the 
shell  coordinates  and  through  the  thickness.  This  generally  is  done  through 
extrapolation  from  previous  solutions. 

(3)  A  subroutine  is  called  within  which,  for  each  of  the  material 
components,  the  stress  corresponding  to  the  assumed  strains  is  determined. 
The  total  stress  for  the  composite  is  then  found. 

J 

°total  Vi  ai 

i=l 

where  J  is  the  number  of  components,  v.  is  the  relative  volume  occupied 
by  component  number  i  v^  =  1.  oj  ,  and  a ^  is  the  yield  stress  for 

component  number  i.  ^ 

(4)  Once  total  strains  and  stresses  are  known,  the  plastic  part 
of  the  strain  increment  can  be  determined  and  added  as  a  pseudo-load 

in  an  elastic  analysis. 

(5)  New  strains  are  computed  and  used  as  estimates.  The  pro¬ 
cedure  is  continued  until  the  corriputer  strains  agree  to  within  a  given  margin 
with  the  estimated  strains. 

The  following  operations  are  performed  in  the  above  referenced  subroutine: 

(1)  Information  about  material  properties  is  transferred  into 

the  subroutine  together  with  the  estimated  strain  increments  (Ae^  ,  A*^  ,  and 
Ay)  and  stresses  at  the  end  of  the  previous  load  step  (o^  ,  ,  y). 

(2)  New  stresses  are  compi  ted  under  the  assumption  that  the 
load  step  is  elastic. 

CT1  =  °i  +— ^2  (A61+  vAe2) 

1  -  v 

a2  ~  °2  +  — ^Ae2  +  (8<?) 

1  -  v 

t  -  t  f  EAy  /  [2  {1  +  vj  J 
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where  k  is  the  ellipse  ratio  for  the  assumed  yield  surface  (usually,  k  =  3). 

2 

(4)  If  Orn  is  less  than  a ^  ,  the  load  step  is  elastic  in  this  com¬ 
ponent  (loading  or  unloading).  If  this  is  the  case  for  all  components,  there 
are  no  psuedo  loads  caused  by  plastic  strain  increments  and  the  calculations 
for  the  load  step  are  concluded. 

2 

(5)  If  is  larger  than  for  some  component,  the  step  is 
at  least  partly  inelastic  for  this  component.  As  we  have  assumed  ideal 
plasticity  the  stresses  can  be  computed  from  the  conditions  that 

2  ,  2  x  ,2  2  2  ,n1, 

al  +  a2  "  ala2  +  k  T  =  aY  (91) 


where 


=  °i  + 


-  Aej3  +  v  (A«2  -  Ae^)J 


°2  a2  + 


-  Ae^  +  v  (A«1  - 


T  =  T  1  •  4vp] 


and  that  the  plastic  flow  is  perpendicular  to  the  yield  surface 


2ql  -  °2 
2a2  -  aL 


2ql  ~  °2 
2k2  7 


After  the  stresses  have  been  determined  in  the  components,  the  avenge 
stress  in  the  composite  is  found  reauily.  As  the  elastic  constants  are 
the  same  for  all  components,  the  plastic  part  of  the  strain  increment 
(i.  e.  ,  the  pseudo  loads),  can  easily  be  obtained.. 
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The  acquired  ability  to  handle  cases  with  inelastic  behavior 
is  demonstrated  in  one  of  the  examples  discussed  in  the  following  section. 


Section  5.  0 

RESULTS  OF  SAMPLE  CASES 


In  the  follov'ing  section  are  presented  some  numerical  results 
obtained  through  exercise  of  the  program.  The  examples  were  chosen 
such  that  the  recent  additions  to  the  program  could  be  verified. 

5.1  Cylinder  With  Rectangular  Cutout 

Analytical  and  experimental  results  for  cylindrical  shells  with 
rectangular  cutouts  were  reported  earlier  in  Reference  3.  The  benefit 
derived  from  the  use  of  a  variable  mesh  spacing  has  been  evaluated  by  re¬ 
examining  this  cylinder  problem.  The  cylinder  has  two  diametrically  oppo¬ 
site  cutouts  and  a  radius -to-thickness  ratio  of  400.  It  was  reported  in 
Reference  3  that'a  reasonably  accurate  analysis  for  such  cylinders  would 
require  excessive  computer  time.  Numerical  results  for  a  uniform  finite 
difference  net  v-ith  9  points  in  the  axial  and  20  points  (9  X  20)  in  the  circum¬ 
ferential  direction  (also  presented  in  Ref.  3)  are  shown  here  in  Figure  6. 
Due  to  improvements  in  the  efficiency  of  the  computer  program,  it  is  now 
possible  to  obtain  much  better  results  even  with  constant  grid  spacing. 
Curve  B  is  obtained  with  a  finer  net  (16  x20),  A  finite  difference  mesh 
was  designed  also  in  which  the  minimum  grid  spacing  is  identical  to  ■ 
that  used  for  Curve  B,  but  which  gradually  increases  away  from  the  cut¬ 
out  until  it  is  approximately  doubled.  The  displacements  corresponding 
to  this  analysis  are  practically  identical  to  those  obtained  by  use  of  grid 
with  constant  spacing,  but  'he  computer  time  is  reduced  by  about  40%. 

Curve  C  was  determined  by  use  of  a  mmimum  grid  spacing  of  0.  2 
inch  at  the  edge  of  the  cutout.  Moving  away  from  the  cutout  the  spacing 
increases  by  a  factor  of  1.  2  from  one  mesh  to  the  next  until  the  maximum 
grid  size,  of  0.  6  inch  is  obtained.  For  Curve  D  the  minimum  si".e  is  0.  1? 
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Figure  6  Load  Displacement  Curves  for  a  Circular  Cylindrical  Shell  with  Cutouts 


inch,  the  factor  is  1.  5  and  the  maximum  size  is  again  0.  6  inch.  The 
results  obtained  by  use  of  the  latter  mesh  appear  to  be  in  very  good 
agreement  with  the  experimental  results.  The  computer  time  correspond¬ 
ing  to  the  determination  of  one  of  these  curves  is  approximately  0.  5  hours 
(UNIVAC  1108).  For  analyses  with  even  finer  mesh  sizes,  therefore,  the 
analysis  was  restricted  to  loads  below  845  pounds.  The  results  in  Table  I 
show  that  additional  refinement  of  the  mesh  would  not  substantially  change 
the  results  shown  in  Curve  D. 


Table  I 


Displac  ement 

w^  at 

P  =  845  lbs 

Net 

Min.  Spacing 

Factor 

Max.  Spacing 

W1 

D  (13  x  21) 

.  12 

1.  5 

.  60 

. 00877 

E  (18  X25) 

.  12 

1.  2 

.60 

.00850 

F  (21  x35) 

.12 

1.2 

.  30 

.  00858 

G  (21  x28) 

.  08 

1.  2 

.  60 

.00873 

5.  2  The  Pincned  Cylinder 

The  case  of  a  pinched  cylinder,  Figure  7,  was  also  analyzed  to 
further  demonstrate  the  advantages  of  the  variable  grid  capability.  Lateral 
displacements  computed  from  a  linear  analysis  are  shown  versus  the  cir¬ 
cumferential  coordinate  in  Figure  8  and  versus  the  axial  coordinate  in 
Figure  9.  The  curves  are  for  a  converged  solution,  corresponding  to  a 
variable  spacing  grid  with  17  points  in  the  axial  and  26  points  in  the  circum¬ 
ferential  directions  (17  X  26).  These  results  are  in  good  agreement  with 
results  for  the  same  case  shown  in  Reference  16.  Discrete  values  of  the 
solution  are  given  for  the  two  coarser  nets  (A  and  B)  which  are  shown  in 
Figure  10.  It  can  be  seen  that  the  use  of  the  net  with  variable  spacing, 

Grid  A,  leads  to  results  which  are  at  least  of  the  same  quality  as  those 
obtained  with  the  uniform  net,  Grid  B.  The  computer  time  corresponding 
to  the  analysis  with  Grid  B  is  approximately  five  times  the  time  for 
analysis  with  Grid  A. 
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Figure  7  Pinched  Cylinder 


* 

Figure  8  Lateral  Displace¬ 
ment  for  Pinched 
Cylinder 


Figure  9  Lateral  Displacements 
for  Pinched  Cylinder 
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GRID  A  (8x12) 


GRID  B  (12x18) 


Figure  10  Finite.  Difference  Grids 
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5.  3  Cylinders  With  a  Circular  Cutout 


A  circular  cylinder  was  analyzed  for  collapse  under  uniform  end 
shortening.  Its  geometrical  properties  were:  length  9  inches, 
radius  6  inches,  thickness  0.  06  inch.  At  its  midlength,  it  had  two  dia¬ 
metrically  opposite  circular  cutouts,  each  of  radius  2.  35  inches.  Young*  s 

7 

modulus  was  set  to  10  psi  and  Poisson's  ratio  to  0.  3.  Due  to  symmetry, 
only  half  the  length  and  one  quarter  of  the  circumference  of  the  shell  was 
considered.  A  finite  difference  net  was  chosen  with  15  axial  and  19  circum¬ 
ferential  stations  (15  x  19).  The  net  is  shown  in  Figure  11.  The  analysis  in¬ 
dicates  collapse  (a  maximum  load)  for  an  end  shortening  of  .  0209  inch.  The 
load  maximum  is  16,  740  lbs  or  66,  960  lbs  for  the  complete  cylinder. 

The  difference  between  the  displacements  for  two  adjacent  solutions 
close  to  the  point  of  collapse  represents  the  collapse  mode  for  the  case. 

In  Figure  12c  is  shown  how  these  incremental  displacements  vary  with  the 
angular  coordinate  (see  Fig.  11).  Figure  12b  shows  the  lateral  displace¬ 
ment  increments  at  the  meridian  8  =  57°.  The  displacements  at  the  edge 
of  :he  cutout  (8  =  22.  5°)  and  at  0  =  57°  are  shown  as  functions  of  the  axial 
load  in  Figure  13.  While  the  largest  displacement  occurs  at  the  cutout  edge, 
the  displacement  at  9  -  57°  is  growing  faster  indicating  that  "buckling" 
occurs  away  from  the  cutout  where  the  axial  stresses  are  higher. 

5.4  Shells  with  Elliptic  Cross-Section 

For  an  elliptic  cone  the  geometric  constants  occurring  in  the  kine¬ 
matic  relations  are  given  as  an  example  in  Section  2.8  .  These  expressions 
were  used  here  in  an  analysis  of  an  elliptic  cone  with  the  dimensions  shown 
in  Figure  14. 

Numerical  results  were  first  obtained  for  the  special  case  of  an 
elliptic  cylinder  with  a  length  of  1.  0  in.  ,  a  thickness  of  0.  0144  in.  ,  and 

7 

semi-axes  of  1.  75  in.  and  1.  0  in.  (sec  Fig..  14).  Young's  modulus  was  10 
psi  and  Poisson's  ratio  was  0.;3.  The.  cylinder  was  subjected  to  a  uniform 
end  shortening  with  the  edges  free  to  rotate  but  restrained  from  moving  in 
the  radial  and  circumferemial  directions. 
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DISPLACEMENT  AT  MIDLENGTH  (In chit) 

Figure  13  Load  Displacement  Curves  for  Cylinder  with  Circular  Cutouts 


Since  the  **buckling  pattern"  was  expected  to  be  confined  to  the 
ar«>as  of  least  curvature,  it  appeared  that  anti  symmetric  behavior  with 
respect  to  the  normal  plane  through  ®  =  0  (Fig.  14)  could  be  excluded. 

Hence,  the  analysis  was  restricted  to  a  180°  panel  with  symmetry  condi¬ 
tions  enforced  at  g  =  0  ,  tt  .  A  uniform  finite  difference  grid  was  chosen 
with  11  points  in  the  axial  and  29  points  in  the  circumferential  directions. 
Results  obtained  with  finer  grids  indicated  that  use  of  the  chosen  grid  led 
to  accurate  computations  of  the  collapse  load. 

Due  to  the  symmetry  of  the  prebuckling  deformation  about  the  plane 
at  midlength  and  about  the  normal  plane  through  cp  =  tt/2  ,  it  was  necessary 
to  excite  nonsymmetric  deformations  by  the  use  of  small  antisymmetric 
imperfections.  Despite  the  presence  of  these  imperfections,  a  deformation 
pattern  developed  at  collapse  which  was  symmetric  about  both  of  these 
planes.  Therefore,  the  continued  analysis  was  restricted  to  panels  cover¬ 
ing  half  the  cylinder  length  and  one  quarter  of  the  circumference. 

For  the  particular  cylinder  considered  (aspect  ratio  of  1.  75),  it  is 
possible  to  determine  the  critical  load  without  the  use  of  symmetric  (with 
respect  to  the  geometric  symmetry  planes)  imperfections.  As  the  load 
is  increased,  a  very  sharp  maximum  is  found  in  the  load  displacement 
curve  (Figure  14).  Beyond  this  maximum  convergence  cannot  be  obtained, 
hence  the  post-buckling  curve  cannot  be  directly  determined. 

For  an  imperfect  shell,  the  displacement  mode  which  developed 
at  collapse  for  a  perfect  shell  was  used  as  a  guide  in  the  choice  of  a  suitable 
initial  imperfection  mode 

w  —  -  wq  sin  cos  (6  8) 

Load  displacement  curves  were  computed  for  several  different  values  of 
the  imperfection  amplitude  wq  The,  results  are  shown  in  Figure  14. 

The  normal  displacement  at  cp  =  tt/ 2  ,  x  =  L/2  is  shown  as  a  function  of 
the  axial  load  in  Figure  15.  From  Figure  14  it  can  be  seen  that  for  a 
sufficiently  large,  imperfection  amplituoe,  the  first  snarp  maximum  does 


-  54  - 


ylinder  s 


200*0 


VM» 


not  exist- -the  curve  is  smooth  and  it  is  possible  to  find  equilibrium  con¬ 
figurations  in  the  post-buckling  range.  After  such  configurations  have 
been  found,  they  can  be  used  as  starting  values  for  an  analysis  in  which 
the  imperfection  amplitude  is  gradually  changed  until  a  point  is  found  on 
the  post-buckling  curve  for  perfect  shells.  After  such  a  point  is  found,  it  is 
easy  to  establish  post-buckling  load  displacement  curves  for  perfect  shells 
(Figure  14). 

After  the  first  sharp  maximum  the  postbuckling  curve  exhibits  two 
additional  limit  points  which  correspond  to  secondary  buckling.  The  curve 
was  not  pursued  beyond  the  third  maximum  because  the  deformations  are 
then  so  large  that  the  applicability  of  the  basic  equations  is  questionable. 
Also  the  buckle  pattern  is  close  to  the  point  of  maximum  curvature  and 
bifurcation  into  an  antisymmetric  mode  is  likely.  The  normal  displace¬ 
ment  as  a  function  of  the  circumferential  arclength  at  x  =  L/2  is  shown 
in  Figure  16.  Curves  are  given  for  each  of  the  three  limit  points  (A,  B, 

C  on  Fig.  14). 

In  the  neighborhood  of  a  limit  point  the  developing  collapse  or 
buckle  mode  can  be  obtained  as  the  difference  between  displacements  for 
two  neighboring  solutions.  Such  collapse  modes  corresponding  to  each  of 
the  three  points  of  maximum  are  shown  in  Figure  17.  It  can  be  seen  that 
the  point  of  maximum  deflection  in  these  patterns  moves  towards  the 
point  of  maximum  curvature  as  the  end  shortening  increases.  While  the 
primary  buckling  load  is  rather  sensitive  to  imperfections,  it  appears  that 
the  second  maximum  is  not  imperfection  sensitive;  hence,  it  may  be  suit¬ 
able  as  a  design  limit.  Results  similar  to  these  have  been  presented  by 
Kempner,  et  al.  ,  for  oval  shells  (Refs.  17,  18).  However,  Kempner's 
shells  are  not  elliptic  and  a  direct  comparison  is  not  possible. 

A  series  of  elliptic  com  s  was  also  analyzed.  Like  the  cylinders, 
the  cones  were  loaded  through  uniform  axial  shortening.  At  the  two  ends 
rotation  was  unrestrained  but  the  cross  section  was  not  allowed  to  deform. 
Four  different  cases  (including  a  circular  cone)  were  i  lalyzed.  The  aspect 
ratio  of  the  elhpfic  cross  section  was  varied  but  the  semi-axes  of  the 
ellipse  were  chosen  such  that  the  circumference  was  the  same  in  all  cases.' 
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Young's  modulus  was  chosen  to  be  10  psi  and  Poisson's  ratio  was  0.  3. 
All  the  cones  had  the  dimensions  (see  Fig.  18)  t  =  0. 16  in. ,  c  =  16  in. , 
and  d  =  16  in.  The  dimensions  of  the  ellipse  are  shown  in  Table  II. 


Table  II 


Case 

(in..). 

>  b 
(in. ) 

1 

10.  65 

10.  65 

2 

11.  9 

9.  5 

3 

12.  2 

8.  7 

4 

13.  0 

7.4 

The  results  for  the  elliptic  cylinders  indicate  that  an  imperfection 
with  an  amplitude  of  about  one  percent  of  the  shell  thickness1  will  not  signi¬ 
ficantly  change  the  critical  load.  However,  if  this  imperfection  is  included, 
a  less  severe  convergence  criterion  may  be  used.  Consequently,  for 
economy  in  the  analysis  such  an  imperfection  was  included  here.  Figure  19 
shows  how  the  critical  load  varies  with  the  ellipse  ratio  for  elliptic  cones 
of  equal  weight.  The  decrease  in  buckling  load  with  increasing  aspect 
ratio  is  less  drastic  than  is  indicated  by  the  bifurcation  buckling  analysis 
for  oi  al  cylinders  (Ref.  16).  For  the  circular  cone  the  bifurcation  point 
and  the  maximum  coincide  but  for  higher  values  of  the  aspect  ratio  the 
critical  load  is  above  the  bifurcation  point.  The  buckling  mode  for  Case'  3 
(a/b  =  1.  4)  is  shown  in  Figure  20.  Similar  results  were  obtained  in  a  bifur¬ 
cation  buckling  analysis  for  oval  cylinders  by  Kempner,  et  al  (Ref.  16)  . 

It  must  be  emphasized  that  for  all  tl.«  cases  investigated  here  a 
uniform  end  shortening  was  applied  tc  the  shell.  Had  a  uniformly  distri¬ 
buted  axial  load  been  applied  at  both  edges,  the  possibilities  for  redistri¬ 
bution  of  stresses  would  have  been  limited  and  the  performance  of  the  elliptic 
shells  would  have  compared  less  favorably  to  shells  with  circular  cross- 
section. 
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5.  5  A  Pear-Shaped  Cylinder 


In  Figure  21  is  shown  a  cylinder  whose  cross  section  consists  of 
circular  arcs  joined  by  straight  lines.  The  behavior  of  this  shell  subjected 
to  uniform  end  shortening  was  investigated  with  use  of  the  STAGS  code. 

As  this  type  of  shell  is  not  included  among  the  standard  geometries, 
a  subroutine  must  first  be  written  for  computation  of  the  geometrical  con¬ 
stants.  The  general  procedure  recommended  in  Reference  8  for  computation 
of  the  geometrical  coefficients  can  be  greatly  simplified  in  a  case  like  this. 

If  the  arclength  and  che  axial  distance  are  chosen  as  surface  coordinates, 
clearly  the  Lame  coefficients  are  A  =  1,  B  =  1  and  C  =  0.  Also  the  local  radii 
of  curvature  are  directly  given. 

As  seen  from  Figure  22  the  linear  range  in  this  case  represents 
less  than  1/30  of  the  total  load  history  of  the  shell.  The  rapid  change  in 
slope  of  the  load-deflection  curves  at  about  P  =  100  lbs  reflects  the  growth 
in  normal  deflection  (buckling)  of  the  flat  portions  of  the  shell.  Associated 
with  this  growth  in  w  is  a  redistribution  of  the  axial  stress  so  that  the 
curved  segments  begin  to  take  up  a  larger  portion  of  the  total  axial  load  P  . 

As  more  and  more  of  the  axial  load  is  borne  by  the  curved  segments,  che  slope 
of  the  load-end-shortening  curve  increases  until  just  before  collapse,  at  which 
load  the  entire  structure  fails.-  Figures  23  and  24  show  the  circumferential 
distributions  of  normal  outward  displacement  w  and  axial  compression/length 
N  at  the  shell  midlength  for  P  =  2328  lbs.  At  this  load,  both  w  and  N 
are  growing  very  rapidly  with  P  in  the  curved  portions  45  £  9  £  90  and 
-67.  5°  £  9  <;  0°. 

The  rather  complex  behavior  in  this  case  indicates  the  need  for  a 
flexible  strategy  for  calculation  of  collapse  loads  of  shells.  Small  load 
steps  and  frequent  refactoring  of  the  equation  system  matrix  are  required 
in  the  load  region  between  100  and  200  lbs  even  though  the  displacements 
are  relatively  small  in  this  range.  Farther  out  on  the  load-end-shortening 
curve,  where  the  displacements  are  larger,-  rather  large  load  steps  can  be 
used  and  few  refactorings  are  necessary.  Efficient  use  of  the  STAGS  code, 
or  any  code  for  predicting  nonlinear  behavior  of  shells,  requires  a  sophisti¬ 
cated  iteration  strategy  built  into  it  and  a  well-trained  user  to  take  advantage 
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of  this  strategy. 

A  finite  difference  grid  was  used  with  45  circumferential  nodes  and 
9  axial  nodes  covering  1/2  of  the  circumference  and  1/2  of  the  length.  A 
variable  spacing  was  used  such  that  the  gridlines  would  follow  the  inter¬ 
sections  be'ween  flat  and  curved  shell  segments. 

5.  6  Bending  of  Cylindrical  Panels  Under  Point  Loading 

The  STAGS  code  was  applied  in  an  analysis  of  the  behavior  of 
shallow  cylindrical  panels  as  shown  in  Figure25.  The  panels  were  subjected 
to  bending  through  application  of  a  point  load  at  the  midpoint  of  a  panel  sup¬ 
ported  at  the  curved  edges  and  with  the  straight  edges  free.  The  behavior 
of  such  shells  is  expected  to  be  highly  nonlinear.  If  the  load  is  applied 
towards  the  center  of  the  circular  arc,  the  cross-section  will  be  more 
and  more  shallow  with  application  of  load  and  the  result  is  similar  to  the 
well  known  Brazier  effect.  If  the  load  is  directed  away  from  the  center, 
the  free  edges  will  be  under  axial  compression  and  the  shell  will  collapse 
under  a  moderate  load. 

Three  cases  were  considered:  one  with  clamped  edges  loaded  to¬ 
wards  the  center  and  two  with  simply  supported  edges;  one  loaded  towards 
and  one  loaded  away  from  the  center.  Ten  axial  and  nine  circumferential 
stations  were  used.  The  results,  in  terms  of  load  displacement  curves, 
are  shown  in  Figure  26  for  the  shells  with  load  toward  the  center,  and  in 
Figure  27  for  the  shell  loaded  away  from  the  center.  In  the  case  with 
clamped  edges,  collapse  is  prevented  by  the  development  of  axial  mem¬ 
brane  tension.  Collapse  in  the  case  of  simple  support  is  indicated  by 
a  maximum  in  the  load  deflection  curve.  In  a  case  like  this,  i.  e.  ,  when 
the  load  is  stepwise  increased  rather  than  a  displacement,  points  on  the 
curve  cannot  be  computed  through  the  maximum.  At  the  point  of  maximum 
the  equilibrium  configuration  is  unstable  and  hence  the  coefficient  matrix 
has  a  zero  determinant.  This  determinant,  as  obtained  when  refactoring 
was  required,  is  plotted  versus  the  load  in  Figure  28.  In  this  case,  it 
is  much  easier  to  read  the  critical  load  from  the  determinant  plot. 
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5.7 


Inelastic  Buckling  of  Plate 


A  flat  plate  was  considered  which  was  simply  supported  on  two 
opposite  edges  and,  on  the  other  two  edges,  in-plane  displacements  were 
allowed  but  lateral  displacements  and  rotation  were  suppressed.  Axial 
compression  was  introduced  in  the  plate  at  the  simply  supported  edges. 

Plate  dimensions  and  boundary  conditions  are  shown  in  Figure  29. 

In  the  elastic  case  the  bifurcation  buckling  theory  would  be  appli¬ 
cable  and  the  value  of  the  critical  load  can  be  obtained  by  use  of  simpler 
means  than  a  nonlinear  analysis.  However,  application  of  STAGS  also 
gives  information  about  the  plate  behavior  in  the  post-buckling  range. 

The  critical  load  for  the  plate  can  be  established  by  use  of  the  nonlinear 
analysis  if  lateral  displacements  are  triggered  by  small  initial  imperfec¬ 
tions.  As  the  lateral  displacements  grow  very  rapidly  and  if  the  imperfec¬ 
tion  amplitude  is  sufficiently  small,  the  load-displacement  curve  has  a 
sharp  knee  at  the  buckling  load  as  it  is  traditionally  defined.  However,  the 
smaller  the  imperfection  is,  the  sharper  must  the  convergence  criterion  be, 
and  the  more  expensive  is  the  analysis.  For  very  small  imperfections,  it 
would  be  necessary  to  use  double  precision  arithmetic.  Therefore,  advan¬ 
tage  was  taken  of  the  fact  that  buckling  is  followed  by  redistribution  of 
stresses.  The  curve  corresponding  to  the  difference  between  axial  stress 
at  the  edge  and  axial  stress  at  the  center  of  the  plate  has  a  much  sharper 
knee  than  the  load  displacement  curve  has  and  it  is  possible  to  determine 
the  buckling  load  with  larger  values  of  the  imperfections. 

The  method  was  demonstrated  first  for  a  plate  which  was  assumed 
to  remain  elastic  for  any  stress.  A  grid  was  used  with  8  nodes  along 
simply  supported  sides  and  6  nodes  along  the  clamped  sides.  The  ini¬ 
tial  imperfection  was  given  by 

,„-5  TTX  T7V 

w  10  sm  -r —  cos  -=~ 
o  LB 

The  results  for  the  elastic  plate  are  shown  in  Figure  30.  The  plot 

2 

of  u  versus  a  indicates  a  value  of  a  critical  load  of  2800  kg/ cm  whieh 
is  m  close  agreement  with  the  result  from  the-  classical  buckling  analysis 
tor  plates. 
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Figure  28  Determinant  Versus  Load  Factor 
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Figure  29  Plate  in  Compression 
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Figure  JO  Load-Dis 


The  same  analysis  was  also  carried  out  for  a  10  percent  thicker 

plate,  both  elastic  and  inelastic  with  the  stress-strain  curve  for  uniaxial 

compression  represented  by  the  polygon  shown  in  Figure  31.  The  results 

are  shown  in  Figure  32  in  the  form  of  a  load-displacement  curve.  A  plot 

of  the  square  of  the  displacement  gives  a  clearer  indication  of  the  critical 

2 

load.  The  critical  stress  is  found  to  be  2270  kg/cm  corresponding  to  an 
axial  load  of  1100  kg.  The  kink  in  the  cm  ve  for  lateral  displacement  is 
presumably  due  to  the  fact  that  as  the  cori. er  on  the  load  displacement  cur/e 
is  reac  ied  by  the  average  stress,  the  bending  stiffness  drops  with  the  reduc¬ 
tion  in  tangent  modulus. 

As  the  form  of  the  load  displacement  curve  above  the  second  corner 
is  in  this  case  irrelevant,  more  accurate  results  would  have  been  obtained 
if  corner  points  had  been  concentrated  in  the  neighborhood  of  the  critical 
stress. 
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Section  6.  0 

CONCLUSIONS  AND  RECOMMENDATIONS 


Certain  improvements  or  extensions  of  the  STAGS  computer  pro¬ 
gram  are  reported  here.  It  appears  from  solutions  of  several  sample 
problems  that  with  these  extensions  the  STAGS  computer  program  has 
become  a  powerful  tool  fer  the  analysis  of  the  nonlinear  behavior  of  shells 
of  general  shape.  The  use  of  the  energy  method  with  finite  differences 
appears  to  be  attractive.  For  most  shells,  one  of  the  standard  geometry 
routines  can  be  used,  in  which  case  determination  of  the  input  data  general¬ 
ly  is  a  matter  of  only  a  few  minutes.  Through  comparison  with  other  pro¬ 
grams  (Ref.  19),  it  has  been  found  that  the  program  is  efficient  wiih  re¬ 
spect  to  computer  time  and  to  numerical  stability.  Likewise,  the  modi¬ 
fied  Newton-Raphson  method  appears  to  be  the  best  choice  for  the  solu¬ 
tion  of  the  nonlinear  equation  system.  It  has  been  favorably  compared 
to  other  numerical  methods  in  Reference  20.  Finally,  through  application 
to  a  large  number  of  practical  cases,,  some  with  previously  known  solu¬ 
tions,  the  validity  of  the  program  has  been  reasonably  well  established  in 
all  its  aspects..  Under  sponsorship  of  the  NASA  Manned  Spacecraft  Center 
in  Houston,  a  series  of  tests  of  cylinders  with  cutouts  has  been  carried  out 
and  results  have  been  compared  to  analytical  resuits  from  STAGS.  The 
agreement  between  test  results  and  analytical  predictions  is  very  good 
(Ref.  21)., 

In  view  of  the  successful  application  of  the  program,  it  appears 
desirable  that  further  extensions  be  made  For  instance,  it  would  enhance 
the  value  of  the  program  if  the  following  items  were  included. 

Improved  ir.p_t  and  output,  particularly  expanded 

p.ot  ■  -a  p  anility  i~v  lucmc  crapniiu.  cispiav  oi  the 


Further  improvement  in  program  efficiency. 

Input  diagnostic. 

Pre-  and  post-processors  of  data  files. 

Inclusion  of  some  finite  elements  in  the  program 
such  as  bars  and  beams,  which  cannot  be  properly 
represented  in  the  present  program.  Such  a  hybrid 
program  would  combine  the  efficiency  of  the  finite 
difference  analysis  with  the  versa* ility  of  the  finite 
element  analysis. 
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